/*
    KMC Simulation
    Copyright (c) 2012 Filip Rozbořil

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
    
    Based on the original implementation by Matsumoto & Nishimura
    
    Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
	 All rights reserved.
*/

#include "MersenneTwister.h"

void MersenneTwister :: generateState() // generate new state vector
{ 
	for (int i = 0; i < (n - m); ++i)
  	{
		state[i] = state[i + m] ^ twiddle(state[i], state[i + 1]);
	}
	for (int i = n - m; i < (n - 1); ++i)
	{
	    state[i] = state[i + m - n] ^ twiddle(state[i], state[i + 1]);
	}
	state[n - 1] = state[m - 1] ^ twiddle(state[n - 1], state[0]);
	p = 0; // reset position
}

void MersenneTwister :: seed(uint32_t s) 
{  // init by 32 bit seed
	state[0] = s;
	for (int i = 1; i < n; ++i) 
	{
		state[i] = UINT32_C(1812433253) * (state[i - 1] ^ (state[i - 1] >> 30)) + i;
		// see Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier
		// in the previous versions, MSBs of the seed affect only MSBs of the array state
		// 2002/01/09 modified by Makoto Matsumoto
	}
	p = n; // force gen_state() to be called for next random number
}

void MersenneTwister :: seed(const uint32_t* array, int size) 
{ // init by array
	seed(UINT32_C(19650218));
	int i = 1, j = 0;
	for (int kk = ((n > size) ? n : size); kk; --kk) 
	{
    	state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * UINT32_C(1664525)))
	      + array[j] + j; // non linear
	   ++j; j %= size;
	   if ((++i) == n) { state[0] = state[n - 1]; i = 1; }
  	}
  	for (int k = n - 1; k; --k) 
  	{
  		state[i] = (state[i] ^ ((state[i - 1] ^ (state[i - 1] >> 30)) * UINT32_C(1566083941))) - i;
	   if ((++i) == n) { state[0] = state[n - 1]; i = 1; }
  	}
  	state[0] = UINT32_C(0x80000000); // MSB is 1; assuring non-zero initial array
  	p = n; // force gen_state() to be called for next random number
}